\documentclass[12pt,titlepage]{article}
\usepackage[spanish]{babel}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{caratulaMetNum}
\usepackage{float}
\usepackage{subfigure}
\usepackage{listings}
\lstset{language=C,basicstyle=\small\tt,keywordstyle=\bf,tabsize=3,breaklines=true,linewidth=16cm,postbreak={\mbox{$\rightsquigarrow$}},prebreak={\mbox{$\rightsquigarrow$}}}

 \usepackage{a4wide}
% \usepackage{amssymb}
% \usepackage{amsmath}
 \usepackage{enumerate}
% \usepackage[utf8]{inputenc}
% \usepackage[spanish]{babel}
% \parindent = 0 pt
 \parskip = 11 pt
% \usepackage[width=15.5cm, left=3cm, top=2.5cm, height= 24.5cm]{geometry}

\usepackage{color}
\usepackage{url}
\definecolor{lnk}{rgb}{0,0,0.4}
\usepackage[colorlinks=true,linkcolor=lnk,citecolor=blue,urlcolor=blue]{hyperref}

\newcommand{\func}[2]{\texttt{#1}(#2) :}
\newcommand{\tab}{\hspace*{2em}}
\newcommand{\FOR}{\textbf{for }}
\newcommand{\TO}{\textbf{ to }}
\newcommand{\IF}{\textbf{if }}
\newcommand{\WHILE}{\textbf{while }}
\newcommand{\THEN}{\textbf{then }}
\newcommand{\ELSE}{\textbf{else }}
\newcommand{\RET}{\textbf{return }}
\newcommand{\MOD}{\textbf{ \% }}
\newcommand{\OR}{\textbf{ or }}
\newcommand{\NOT}{\textbf{ not }}
\newcommand{\tOde}[1]{\tab \small{\mathcal{O}($#1$)}}
\newcommand{\Ode}[1]{\ensuremath{\small{\mathcal{O}\left(#1\right)}}}
\newcommand{\VSP}{\vspace*{3em}}
\newcommand{\Pa}{\vspace{5mm}}
\newenvironment{pseudo}{\noindent\begin{tabular}{ll}}{\end{tabular}\VSP}

\newenvironment{while}{\WHILE \\ \setlength{\leftmargin}{0em} }{}

\newcommand{\iif}{\Leftrightarrow}
\newcommand{\gra}[1]{\noindent\includegraphics[scale=.70]{#1}\\}
\newcommand{\gras}[2]{\noindent\includegraphics[scale=#2]{#1}\\}
\newcommand{\gram}[1]{\noindent\includegraphics[scale=.50]{#1}}
\newcommand{\dirmail}[1]{\normalsize{\texttt{#1}}}
\newenvironment{usection}[1]{\newpage\begin{section}*{#1}	\addcontentsline{toc}{section}{#1}}{\end{section}}
\newenvironment{usubsection}[1]{\begin{subsection}*{#1}	\addcontentsline{toc}{subsection}{#1}}{\end{subsection}}

\newcommand{\superref}[1]{\textsuperscript{\ref{#1}}}

%\title{{\sc\normalsize Métodos Numéricos}\\{\bf Trabajo Práctico Nº1}}
%\author{\begin{tabular}{lcr}Pablo Herrero & LU & \dirmail{pablodherrero@gmail.com}\\Thomas Fischer & 489/08 & \dirmail{tfischer@dc.uba.ar}\\Kevin Allekotte & 490/08 & \dirmail{kevinalle@gmail.com} \end{tabular}}
%\date{\VSP \normalsize{Abril 2010}}
\begin{document}

\materia{Métodos Numéricos}
\titulo{Trabajo Práctico Nº2}
\subtitulo{Chocan los Planetas...}
%\grupo{Grupo x}
\integrante{Pablo Herrero}{332/07}{pablodherrero@gmail.com}
\integrante{Thomas Fischer}{489/08}{tfischer@dc.uba.ar}
\integrante{Kevin Allekotte}{490/08}{kevinalle@gmail.com}

\abstracto{
	El siguiente trabajo se propone la comparación entre distintos
	métodos de simulación de sistemas gravitatorios
	de varios cuerpos usando resolución de sistemas lineales.
	El objetivo final es calcular las coordenadas para tirar
	un proyectil desde la órbita del planeta Neptuno para
	que el mismo impacte	en la Tierra.
}

\palabraClave{Simulación}
\palabraClave{Gravitación}
\palabraClave{Sistemas Lineales}
\palabraClave{Proyectil}

\begin{titlepage}
\maketitle
\end{titlepage}
\tableofcontents
\newpage

	\begin{usection}{Introducción teorica}

		En la mecánica classica o Newtoniana, el problema de lo $n$
		cuerpos es el problema de predecir el movimiento de cuerpos
		celestes en mutua interacción gravitatoria para cualquier
		instante en el futuro o deducir su movimiento en cualquier
		instante del pasado.

		La acción gravitatoria entre los cuerpos depende de las
		posiciones relativas entre ellos, por lo que las ecuaciones de
		interacción gravitatoria quedan definidas en términos de
		ecuaciones diferenciales.

		El problema que surge al plantear el problema con $3 \leq n$ es
		que las ecuaciones el sistema quedan definidas en función de
		demasiadas variables como para aplicar los métodos de integración
		conocidos hoy en día para resolver las ecuaciones.

		Luego tenemos que encontrar otra forma de calcular las
		variables. En nuestro caso elegimos el método de simulación,
		el cuál introduce errores analíticos por discretizar el tiempo y
		de redondeo por trabajar con la aritmética finita de la
		computadora, pero que tiene la ventaja que el sistema a resolver
		es lineal, para lo cuál conocemos un método.

		Para minimizar estos errores implementamos distintos métodos de
		cálculo numérico de aproximación que validamos con sistemas
		familiares como el sistema solar, del cual sabemos por ejemplo
		que las órbitas son elípticas y bastante estables. Si dibujamos
		las órbitas simuladas podemos ver si se comportan de manera
		similar a lo esperado. Además, tenemos datos que se supone son
		bastante precisos descargados de la web oficial de la NASA para
		comparar.

	\end{usection}

	\begin{usection}{Desarrollo}

		En primer lugar el desarrollo consistio en implementar módulos
		para las operaciones con y entre Matrices y Vectores que
		fuesen necesarios para resolver sistemas lineales. Esto
		implicaba también implementar el algoritmo de eliminación
		Gaussiana con pivoteo parcial para resolver el sistema lineal
		del segundo método de resolución. Los módulos fueron
		implementados en el lenguaje \texttt{C++} para garantizar una
		buena performance.
		
		Escribimos el módulo \texttt{Matriz.h} que implementa una matriz
		con las operaciones:
		
		\begin{tabular}{rl}
			\texttt{T} & Transpuesta\\
			\texttt{resolver} & Solución al sistema $Mx=b$\\
			\texttt{triangular} & Triangula la matriz con el método de eliminación gausseana con pivoteo simple\\
			\multicolumn{2}{l}{
				\texttt{Suma, Resta, Multiplicación, Multiplicación por escalar}
			}
		\end{tabular}
		
		Elegimos implementar la triangulación con el método gausseano
		de pivoteo simple porque es una de las formas más simples para
		resolver un sistema de ecuaciones lineal.
		
		Además escribimos la clase \texttt{Vector3} que implementa un vector
		en $\mathbb{R}^3$ con las operaciones básicas como \texttt{Suma, Resta
		Producto Escalar, Producto Vectorial, Normalizar}, etc.

		Luego procedimos a implementar un entorno de simulación
		(\texttt{C++}) y un pequeño script para plotear y poder
		visualizar los resultados (\texttt{Python}). 

		Con todo el backend funcionando, procedimos a implementar
		gradualmente los distintos algoritmos de solución de sistemas
		lineales, intercalando la producción de cada uno, con fases de
		validación que consistían en simular los sistemas propuestos
		en el enunciado, y cuyos resultados se pueden observar en la
		sección \textbf{resultados}.
		
		Para validar los métodos utilizados en la simulación corrimos pruebas
		que consistían en casos simples como por ejemplo Sol-Tierra-Luna, o el
		sistema solar y observamos si las trayectorias de los cuerpos se comportaban
		como esperadas. Por ejemplo comprobamos si la órbita de la tierra alrededor del
		sol era aproximadamente circular y si recorría una vuelta entera en 365 días.
		Ademas vimos si la trayectoria de la luna parecía acompañar a la tierra.
		Observamos que con un $dt$ de una hora los planetas se comportaban suficientemente
		parecido a lo que esperábamos.

		Validados los métodos de simulación, tuvimos que idear un
		algoritmo para calcular el vector de velocidad inicial de
		nuestro proyectil dirigido a la tierra, de tal manera que le pegase.
		La idea fue hacer una especie de búsqueda local, ya que
		suponíamos que de no colisionar con otro planeta en la
		trayectoria hasta la tierra, si pasabamos con una trayectoría
		mas cerca de la tierra que con otra, podíamos reducir la
		búsqueda de la solución a un espacio localizado cerca de la
		primera. Con un espacio o una trayectoria nos referimos a una
		velocidad inicial, la que en nuestro código esta determinada por
		dos ángulos de disparo. Uno alrededor el eje $\hat{z}$, sobre el
		plano $\hat{x}\hat{y}$ ($\phi$) y otro alrededor del eje $\hat{y}$, sobre
		el plano $\hat{x}\hat{z}$ ($\theta$).
		
		Dadas las condiciones iniciales de un sistema (posiciones y velocidades de todos los cuerpos), podemos simular este sistema con alguno de los métodos y calcular las posiciones para algun tiempo en el futuro.
		Nuestro objetivo es encontrar la velocidad inicial de nuestro proyectil para que esté a menos de $10^{-4}$ AU de la tierra en algun instante futuro.
		Definimos entonces una función (\texttt{mindist}) que dadas las condiciones iniciales, un cuerpo \texttt{proyectil} y un cuerpo \texttt{target}, nos devuelva la distancia mínima a la que se encontrarán estos cuerpos en algun futuro acotado.
		Nos restringimos a los casos en los que la distancia es monótonamente decreciente, y consideramos sólo los tiempos para los cuales los cuerpos se acercan. (Esto es una desicion de implementación, ya que no estaríamos encontrando las soluciones en las que el proyectil se aleja y luego se vuelve a acercar).
		
		Como posición de lanzamiento elegimos un punto aproximadamente sobre la
		órbita de neptuno, esto es, a unos 30AU del sol, sobre el plano $xy$.
		Por ejemplo, los lanzamientos de prueba tenían como posición inicial el punto $(30,0,0)$.
		
		Tenemos entonces que minimizar esta función (o encontrar un valor para el cual sea menor que $10^{-4}$).
		Para esto elegimos una velocidad inicial a mano (aproximadamente hacia el cuerpo \texttt{target}) y corremos un algoritmo que va ajustando la velocidad para minimizar la función \texttt{mindist}. [Más detalles en \textbf{Discusión}]

		\begin{figure}[H]
			\centering
			\includegraphics[scale=0.6]{mindist3.png}
			\caption{Esta es la forma de la función \texttt{mindist} para ciertas condiciones iniciales, donde el dominio es $\phi\times\theta$ (la velocidad inicial del proyectil) y la imagen es la distancia mínima.}
			\label{fig:mindist}
		\end{figure}	

	\end{usection}
	\begin{usection}{Resultados}

	Los experimentos que realizamos fueron los propuestos por el enunciado del trabajo,
	además de algunos gráficos de análisis de error con respecto a los datos de referencia (NASA)
	y simulaciónes de validación de la luna alrededor de la tierra.

	Para cada simulación, se tomaron los datos iniciales de posición y velocidad en la fecha propuesta
	de los datos de la NASA.
	
	\label{sec:res}
		\input{resultados_ej1_m1.tex}
		\input{resultados_ej1_m2.tex}
		\input{resultados_ej1_m3.tex}
		\input{resultados_ej1_error.tex}
		\input{resultados_ej1_relativo.tex}
		\input{resultados_ej2_m1.tex}
		\input{resultados_ej2_m2.tex}
		\input{resultados_misil.tex}
		\input{resultados_bomba.tex}
	\end{usection}
	
	\begin{usection}{Discusión}
		
		
		\begin{usubsection}{Algoritmo de búsqueda del mínimo \texttt{mindist}}
		Para cada trayectoria inicial, definimos una grilla de $3 \times 3$
		puntos, donde el punto central era ella misma, y los otros
		puntos son las rotaciones de la velocidad sobre los ejes
		anteriormente mencionados por dos ángulos, respectivamente, de
		un módulo inicial que se reduce a la mitad en cada iteración de
		la búsqueda, y luego disparámos proyectiles en todas las 9
		direcciónes.
		Así, aunque en algunos casos asintóticamente,
		podemos llegar a cualquier punto del espacio
		contenido en la grilla inicial (la definida en la primera
		iteración) si alli se encuentra una solución.
		A la grilla inicial le ponemos un ancho inicial grosero, de
		manera de incluír casi necesariamente a alguna solución.

		\begin{figure}[H]
			\centering
			\includegraphics[scale=0.5]{img/grid.png}
			\caption{
				En la imagen vemos una representación gráfica de los
				sucesivos espacios de búsqueda local. el orden de las
				iteraciones es rojo $\rightarrow$ azul
				$\rightarrow$ verde.
			}
			\label{fig:grilla}
		\end{figure}		

		Para elegir la trayectoria inicial del algoritmo, calculamos
		primero aproximadamente (simulando y visualizando) cuanto
		tiempo tardaba el proyectil en llegar hasta la órbita terrestre
		desde la órbita de neptuno. Luego calculamos aproximadamente
		la posición de la Tierra en esa cantidad de tiempo después de la
		fecha inicial de lanzamiento, tomando los datos de la posición
		de ese día de la información que provee la NASA.

		También de esos datos, elegimos un punto arbitrario en la
		órbita de Neptuno para el lanzamiento del proyectil. Como
		dirección de la velocidad inicial del proyectil, elegimos el
		vector resultante de restar la posición inicial del proyectil a la
		posición aproximada del impacto sobre la órbita terrestre.
		El módulo de la velocidad inicial, es el propuesto por el
		enunciado del problema.

		Luego de probar para un par de posiciones y velocidades
		iniciales distintas, encontramos soluciones donde pudimos
		destruír la tierra. Las mismas están detalladas en la sección
		siguiente.
		\end{usubsection}
		
		\begin{usubsection}{Gauss}
			Elegimos resolver los sistemas lineales usando el método de Gauss porque es fácil de implementar y es un poco más eficiente que los métodos de factorización.
			No se justifica usar \texttt{LU} porque la ventaja principal es que hay que factorizar la matriz una sola vez para resolver muchos sistemas, pero nosotros resolvemos siempre sistemas con diferentes matrices.
		\end{usubsection}

		\begin{usubsection}{$\Delta t$}
			Un problema que observamos durante las simulaciones, cuando tratábamos de pegarle con el proyectil a la tierra, es que el $\Delta t$ no era lo suficientemente chico para lograr la presicion de impacto requerida.
			De un instante de tiempo al otro, el proyectil estaba adelante y luego atrás de la tierra, siendo la distancia entre los cuerpos nunca suficientemente chica.
			Por esto en la función \texttt{mindist} alteramos el $\Delta t$ dinámicamente a medida que el proyectil se acercaba al \texttt{target}, para tener más presicion cuando se encontraban cerca.
		\end{usubsection}

		\begin{usubsection}{Error}
		
			En los gráficos respectivos a la tierra y la luna en la figura \ref{res_ej1_err}, vemos
			como el error evoluciona de manera pareja en función del tiempo para el método 1.
			Esto parece indicar un buen funcionamiento del algoritmo, ya que no hay saltos bruscos en la gráfica del error.
			
			Los errores se deben seguramente a la diferencia de precisión en el cálculo con respecto a los cómputos de la NASA,
			y se incrementan contínuamente a causa del error acumulado.
			Podemos también ver que el comportamiento del error
			se ve reflejado en los gráficos correspondientes
			de las órbitas de los planetas mostradas anteriormente, y viceversa.
			
			En los gráficos correspondientes al método 2 podemos ver que se producen osilaciones
			que curiosamente tienen un peíodo aproximado de 365 días.
			Si buscamos la causa en las órbitas, podemos decir que a medida que
			la tierra y por ende la luna se alejan de su punto de partida,
			se alejan también de la órbita ``correcta'',
			que suponemos es muy parecida a la calculada por la NASA.
			Cuando pasado la mitad del período de orbitación terrestre,
			los cuerpos comienza a acercarse de nuevo a su punto de partida,
			la órbita simulada se va ajustando de nuevo a la de la NASA.
			O sea que la elipse de la órbita es un poco mas grande que la real,
			pero sin embargo tienen un punto en ``común''.
			Igual vemos que no se ajusta del todo, ya que en cada período
			se incrementa la distancia a la órbita ideal.

			Esto nos lleva a pensar que tal vez halla una falla en la implementación del algoritmo correspondiente al método 2,
			posiblemente al armar la matriz con los datos o en el algoritmo de Gauss para solución de ecuaciones lineales.
			No obstante, vemos que a medida que achicamos el $\Delta t$ el error disminuye considerablemente, por lo cuál tal vez
			el método simplemente no sea del todo adecuado para trabajar con esta granularidad,
			y disminuyendola más (lo cuál requiere mucho tiempo de cómputo, por eso no tenemos resultados)
			se logren resultados buenos.

			El método 3 nos resulto demasiado costoso para calcular
			durante un período de tiempo significativo como para hablar
			de la evolución del error.
			
		\end{usubsection}
		
	\end{usection}
	
	\begin{section}*{Conclusiones}	\addcontentsline{toc}{section}{Conclusiones}

		Llegamos a la conclusión que, dados los métodos implementados,
		el mejor en la práctica para el ejercicio de destruír la tierra, resulto ser el método 1, ya que,
		los otros solo resultaban muy buenos (mejores incluso en cuanto a los resultados finales que el método 1)
		si el $\Delta t$ era significativamente muy bajo,
		lo que para nuestro algoritmo de búsqueda local era muy malo en cuanto a tiempo de ejecución.

		No obstante, más alla de la sensibilidad del método 2 a los $\Delta t$ altos,
		el método de triangulación resulto ser bastante mas rápido de lo que esperabamos, en la práctica.

		Un detalle que notamos fue que la precisión relativa entre el método 1 y 2 de la simulación, a medida que el $\Delta t$ se iba haciendo mas chico,
		era muy similar. Lo que concluímos fue que a pesar de que los 2 algoritmos son bastante distintos,
		y el segundo pareciera ser mucho más sofisticado, las aproximaciones que hacemos en ambos casos son siempre lineales,
		lo que parece justificar nuestras observaciones sobre la similitud de los métodos antes mencionada,
		y nos dice que a partir de un cierto $\Delta t$ suficientemente chico (ya que el método 2 es mas sensible con $\Delta t$ grande) las aproximaciones del segundo método
		no van a ser mejores que las del primero.

		El método 3 parece arreglar este problema, aproximando varias derivadas en un mismo $\Delta t$,
		hasta que la diferencia entre dos derivadas sucesivas se haga mas chica que un umbral.
		Sin embargo esto es muy parecido, y practicamente igual de eficiente que achicar todavía mas el $\Delta t$ para los métodos anteriores.

	\end{section}
	
	\begin{usection}{Apéndices}
		\begin{usubsection}{Apéndice A: Enunciado}
			\input{enunciado.tex}
		\end{usubsection}
		
		\newpage
		\begin{usubsection}{Apéndice B: Codigos Fuente}

\textbf{Algunos defines}
\begin{lstlisting}
#define G 0.000148818071108351462554986428198016585385
#define F(i,j,x) (-G*Cuerpos[i].m*Cuerpos[j].m*(y[3*N+3*i+x]-y[3*N+3*j+x])/(dist*dist*dist))
#define AC(i,j,x) (-G*Cuerpos[j].m*(y[3*N+3*i+x]-y[3*N+3*j+x])/(dist*dist*dist))
#define yX(i) y[3*N+3*(i)]
#define yY(i) y[3*N+3*(i)+1]
#define yZ(i) y[3*N+3*(i)+2]
#define XYZ(i) V3(yX(i),yY(i),yZ(i))
#define VEL(i) V3(y[3*(i)],y[3*(i)+1],y[3*(i)+2])
#define sq(x) ((x)*(x))
\end{lstlisting}

\bigskip
\textbf{Extracto de \texttt{main}, que minimiza \texttt{mindist}}
\begin{lstlisting}
int intento=0;
while(min.first>1e-4 && intento<26){
	intento++;
	for(int ii=-1;ii<=1;ii++) for(int jj=-1;jj<=1;jj++){
		Cuerpos[misil_index].v = pdir.rotate(ii*span,jj*span); // Calculo direccion
		y=makeY(); // Rehago el vector y
		pair<long double,int> new_min = mindist(y,misil_index,target_index);
		if( new_min.first<min.first || ( new_min.first==min.first && new_min.second<min.second ) ){
			min = new_min;
			mindir = Cuerpos[misil_index].v;
		}
	}

	pdir=mindir;
	span*=.6;
	clog << "mindist: " << min.first << " span: " << span << " dir: " << pdir << " it: " << min.second << endl;
}
\end{lstlisting}

\bigskip
\textbf{Función \texttt{mindist}}
\begin{lstlisting}
pair<long double,int> mindist(const Vn& y_in, int obj, int target){
	/* Calcula la minima distancia a la que le pasa el proyectil a la tierra */
	// Devuelve la distancia minima que alcanzan los objetos obj y target mientras se esten acercando y mientras no supere el tiempo de simulacion
	double dtbak=dt;
	Vn y(y_in);
	long double d=(XYZ(obj)-XYZ(target)).norm();
	long double dans;
	int i=0;
	do{
		dans=d;
		y=next(y);
		d=(XYZ(obj)-XYZ(target)).norm();
		if(VEL(obj).norm()*dt>d*.005){ dt*=.5; /*clog << "dt: " << dt << endl;*/}
	}while(d<dans && ++i<resolution);
	dt=dtbak;
	return pair<long double,int>(dans,i);
}
\end{lstlisting}

\bigskip
\textbf{Funciones que calculan la matriz Df}
\begin{lstlisting}
Matriz dFdx(int i, int j, const Vn& y){
	long double d=sqrt( sq(yX(i)-yX(j)) + sq(yY(i)-yY(j)) + sq(yZ(i)-yZ(j)) );
	long double d3=1/(d*d*d);
	long double d5=3/(d*d*d*d*d);
	long double dx=yX(i)-yX(j);
	long double dy=yY(i)-yY(j);
	long double dz=yZ(i)-yZ(j);
	Matriz dFdx(3,3,0);
	dFdx(0,0)=d3-d5*dx*dx; dFdx(0,1) = -d5*dx*dy; dFdx(0,2) = -d5*dx*dz;
	dFdx(1,0) = -d5*dx*dy; dFdx(1,1)=d3-d5*dy*dy; dFdx(1,2) = -d5*dy*dz;
	dFdx(2,0) = -d5*dx*dz; dFdx(2,1) = -d5*dy*dz; dFdx(2,2)=d3-d5*dz*dz;
	return -G*Cuerpos[i].m*Cuerpos[j].m*dFdx;
}

Matriz Df(const Vn& y){
	Matriz res(6*N,6*N,0);
	// lleno A21
	forn(i,3*N) forn(j,3*N) if(i==j) res(3*N+i,j) = 1;
	// lleno A12
	forn(i,N) forn(j,N) {
		Matriz S(3,3,0);
		if(i==j){
			forn(k,N) if(k!=i) S += dFdx(i,k,y);
			S *= (1./Cuerpos[i].m);
		}else{
			S = (1./Cuerpos[i].m)*dFdx(i,j,y);
		}
		forn(ii,3) forn(jj,3) res(3*i+ii,3*N+3*j+jj) = (long double)S(ii,jj);
	}
	return res;
}
\end{lstlisting}

\bigskip
\textbf{Funcion f(y) ($\frac{\delta y}{\delta t}$)}
\begin{lstlisting}
Vn f(const Vn& y){
	Vn res(6*N,0);
	forn(i,N){
		long double sumx=0,sumy=0,sumz=0;
		forn(j,N) if(j!=i){
			long double dist=(V3(y[3*N+3*i],y[3*N+3*i+1],y[3*N+3*i+2])-V3(y[3*N+3*j],y[3*N+3*j+1],y[3*N+3*j+2])).norm();
			sumx+=AC(i,j,0);
			sumy+=AC(i,j,1);
			sumz+=AC(i,j,2);
		}
		res[3*i]  =sumx;
		res[3*i+1]=sumy;
		res[3*i+2]=sumz;
	}
	forn(i,3*N) res[i+3*N]=y[i];
	return res;
}
\end{lstlisting}

\bigskip
\textbf{Método 2 (Taylor)}
\begin{lstlisting}
Matriz Taylor(const Vn& y){
	Matriz Dfy = Df(y);
	Matriz A( Matriz::ID(6*N) - dt*Dfy );
	return A.resolver( y + dt*( f(y) - Dfy*y ) );
}
\end{lstlisting}

\bigskip
\textbf{Método 3 (Iterativo)}
\begin{lstlisting}
Vn MetodoIterativo(const Vn& y, const long double dx){
	Vn w0(y); Vn w1(y);
	long double d0 = INF; long double d1 = INF;
	int i=0;
	do{
		w0 = w1;
		Matriz Dfw = Df(w0);
		Matriz A( Matriz::ID(6*N) - dt*Dfw );
		Matriz B( y + dt*( f(w0) - Dfw*w0 ) );
		w1 = A.resolver( B );

		d0 = d1;
		d1 = dist(w1-w0);
		i++;
	}while( d1<d0 && d1<dx );
	return w0;
}
\end{lstlisting}

\bigskip
\textbf{Resolucion de sistemas lineales (Gauss)}
\begin{lstlisting}
Matriz Matriz::resolver(const Vector& b){
	double U2[n*m]; forn(i,n*m) U2[i] = M[i];
	int P2[n]; forn(i,n) P2[i] = i;

	//double res[n]; forn(i,n) res[i] = b.elem(i,0);
	Vector res(b);

	forn(k,n-1) triangular(res,U2,P2,k);

	fornr(i,n){
		forsn(j,i+1,n) res[i] -= U2[i*m + j]*res[j];
		res[i] /= U2[i*m + i];
	}
	return res;
}

void Matriz::triangular(Vector& res, double* U2, int* P2, const int k){
	int f = k;
	double p = U2[f*m + k];

	// elijo la fila pivote donde el k-esimo elem es maximo
	forsn(i,k,n) if( abs(U2[i*m + k]) > abs(p) ){ f = i; p = U2[f*m + k]; }

	// swapeo las filas f y k en U y res
	forn(j,n){
		double temp_U2 = U2[f*m + j];
		U2[f*m + j] = U2[k*m + j];
		U2[k*m + j] = temp_U2;
		double temp_r = res[f];
		res[f] = res[k];
		res[k] = temp_r;
	}

	// actualizo el vector de perumtacion
	int temp = P2[f];
	P2[f] = P2[k];
	P2[k] = temp;

	// triangulo las filas k+1 -> n
	forsn(i,k+1,n){
		// elijo el m_ij
		double a = U2[i*m+k]/p;
		// la k-esima columna es 0
		U2[i*m+k] = 0;
		// calculo el valor de las columnas k+1 -> n
		forsn(j,k+1,n) U2[i*m+j] -= a*U2[k*m+j];
		res[i] -= a*res[k];
	}
}
\end{lstlisting}

		\end{usubsection}
	\end{usection}
	
	\begin{usection}{Referencias}

	\begin{enumerate}
		\item \label{ref:nasa} \texttt{http://ssd.jpl.nasa.gov/horizons.cgi} \\ Datos de la NASA sobre las posiciones y velocidades de los planetas
		\item \label{ref:matplotlib} \texttt{http://matplotlib.sourceforge.net/contents.html} \\ Referencia de matplotlib, para los gráficos
	\end{enumerate}	
	
	\end{usection}

\end{document}

